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Abstract 

This  paper  introduces  a  new  fanaily  of  second-order  methods  for  solving 
the  index-2  DAE  equations  of  motion  for  flexible  mechanism  dynamics.  These 
methods,  which  extend  the  a-methods  for  ODEs  of  structural  dynamics  to 
DAEs,  possess  numerical  dissipation  that  can  be  controlled  by  the  user.  Con¬ 
vergence  and  stability  analysis  is  given  and  verifies  that  the  DAE  a-methods 
introduce  no  additional  oscillations  and  preserve  the  stability  of  the  underlying 
ODE  system.  Convergence  of  the  Newton  iteration,  which  can  be  a  source  of  dif¬ 
ficulty  in  solving  nonlinear  oscillatory  systems  with  large  stepsizes,  is  achieved 
via  a  coordinate- split  modification  to  the  Newton  iteration.  Numerical  results 
illustrate  the  effectiveness  of  the  new  methods  for  simulation  of  flexible  mech¬ 
anisms. 


1  Introduction 

The  numerical  solution  of  flexible  multibody  systems  is  required  for  geometrically 
nonlinear  dynamic  analysis  of  articulated  structures.  These  are  structures  in  which 

*The  work  was  partially  sponsored  by  the  Army  High  Performance  Computing  Research  Center 
under  the  auspices  of  the  Department  of  Army,  Army  Research  Laboratory  cooperative  agreement 
numbers  DAAH04-95-2-0003/contract  number  DDAH04-95-C-0008,  and  by  ARO  contract  number 
DAAH04-94-6-0208  and  DAAH04-94-6-0213,  the  content  of  which  does  not  necessarily  reflect  the 
position  or  the  policy  of  the  government,  and  no  official  endorsement  should  be  inferred. 
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kinematic  connections  permit  large  relative  motion  between  components  that  undergo 
small  elastic  deformation.  A  source  of  difficulty  in  the  solution  of  flexible  multibody 
equations  of  motion  is  the  coupling  between  the  elastodynamic  equations  and  the 
gross  motion.  Simulation  of  flexible  multibody  systems  has  been  an  active  research 
topic  for  the  last  two  decades.  Many  of  the  developments  in  flexible  multibody  sys¬ 
tems  have  been  implemented  in  the  multibody  dynamic  analysis  codes  [19,  21].  For 
such  systems,  an  important  feature  of  the  solution  is  the  nonlinear  oscillations  in¬ 
duced  by  the  elastodynamics  equations.  Moreover,  since  the  governing  equations  of 
flexible  multibody  systems  axe  often  modeled  using  algebraic  constraints,  the  numeri¬ 
cal  solution  of  differential-algebraic  equations  (DAE)  is  required  [1,  27].  This  paper  is 
concerned  with  the  numerical  solution  of  the  DAE  of  constrained  flexible  multibody 
motion  by  applying  the  family  of  a-methods  for  linear  structural  dynamics  problems. 

Recent  developments  in  the  simulation  of  multibody  systems  have  created  a  variety 
of  methodologies  for  efficient  automation  of  the  process  of  constructing  the  highly 
nonlinear  equations  of  motion  of  constrained  multibody  mechanical  systems  [10], 

q-  +  =  0  (1.1a) 

g(q,t)  =  0  (1.1b) 

where  q  is  the  generalized  coordinates,  A  is  the  Lagrange  multipliers  that  couple  the 
kinematic  constraints  g  with  the  Newton-Euler  equations  of  motion,  and  G  =  |^  is 
the  constraint  Jacobian.  While  many  contributions  to  the  literature  have  focused  on 
the  efficient  generation  of  (1.1)  [10,  14,  20,  21],  in  recent  years  there  has  been  progress 
in  the  numerical  solution  of  differential- algebraic  equations  (DAE)  [1,  9,  8].  Solution 
techniques  for  multibody  systems  have  been  proposed  in  [5,  6,  15,  18,  23,  26].  Two 
basically  different  approaches  have  been  proposed  in  the  literature  for  integration  of 
(1.1).  The  first  involves  direct  application  of  stiffly  stable  numerical  integration,  such 
as  backward  differentiation  formulas  (BDF),  to  various  forms  of  DAEs  derived  from 
(1.1).  The  second  approach  is  based  on  the  solution  of  the  state-space  form  of  (1.1) 
that  is  projected  onto  the  constraint  manifold. 

Time  integration  algorithms  for  solving  structural  dynamics  problems  have  been 
developed  since  the  late  1950’s  [17].  General  requirements  and  the  foundations  of 
these  methods  have  been  well-documented  [11,  12,  28,  3].  Although  their  main  appli¬ 
cation  area  is  to  linear  structural  dynamics,  these  methods  can  be  directly  applied  to 
the  initial  value  problems  of  nonlinear  second-order  ordinary  differential  equations 

(ODE), 

«  =  /(9>9.0-  (1-2) 

Accuracy  and  stability  analysis  hold  for  the  numerical  methods,  provided  the  dis¬ 
cretized  nonlinear  equations  have  been  solved  accurately,  i.e.,  within  a  small  enough 
tolerance.  For  example,  the  HHT-a  method  [11]  for  (1.2)  is  given  by 

Un+i  =  (1 -f  a)/„+i  -  a/n  (1.3a) 
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(1.3b) 

(1.3c) 


?n+i  =  9n  +  A»n  +  A^((r  “  ;3)a„  +  /8a„+i) 
fn+l  =  Kn  +  4((1  -  7)“n  +  70i.+l) 

where  the  subscript  denotes  the  corresponding  function  value  at  the  time,  e.g., 
Qn  =  q{tn),  h  =  tn+i  —  tn  is  the  time  Step,  a  €  [— 1,0],  ^  ,  and  7  =  |  —  a.  It  is 

well-known  that  the  HHT-a  family  is  second-order  accurate  and  A-stable.  Numerical 
dissipation  is  maximum  for  oc  =  —0.3,  and  zero  for  o:  =  0.  Controllable  numerical  dis¬ 
sipation  and  unconditional  stability  are  needed  to  deal  with  the  high-frequency  modes 
which  often  result  from  standard  finite  element  spatial  discretization.  For  nonlinear 
oscillations,  these  properties  are  also  required  in  the  solution  of  flexible  multibody 
systems.  Rather  than  using  ad  hoc  mode-selection  processes,  this  approach  is  more 
desirable  because  the  elimination  of  higher  frequencies  is  controlled  by  selection  of 
the  parameters. 

This  paper  is  concerned  with  the  numerical  solution  of  constrained  flexible  multi¬ 
body  motion  by  the  family  of  ot-methods.  Attempts  have  been  made  to  apply  the 
H HT-a-meihod  (1.3)  directly  to  (1.1)  and  its  index-2  form  [2].  However,  these  meth¬ 
ods  are  plagued  by  oscillations  which  are  largely  unphysical.  Here,  we  extend  the 
Q!-method  to  DAEs  in  a  way  which  does  not  introduce  additional  oscillations.  The 
basic  idea  of  the  proposed  method  is  to  use  the  underlying  ODEs  of  DAEs  and  the 
projection  onto  the  constraint  manifold.  The  main  effort  is  in  developing  numerical 
schemes  which  exploit  the  structure  of  multibody  dynamic  equations,  to  maintain 
the  stability  and  reduce  the  computational  cost  as  compared  to  a  straightforward 
application.  To  treat  the  internal  oscillations  in  the  constrained  systems,  we  use  a 
modified  Newton-type  iteration  based  on  the  Coordinate- Split  formulation  of  multi¬ 
body  systems  [26]. 

In  Section  2,  we  show  the  direct  application  of  the  a-method  to  nonlinear  multi¬ 
body  DAEs.  A  projection  step  to  enforce  the  position  and  velocity  constraints  is 
used  for  maintaining  the  stability  and  accuracy  of  the  numerical  solution.  A  con¬ 
vergence  analysis  of  this  approach  is  given.  In  Section  3,  we  describe  an  effective 
form  of  the  a-method  for  constrained  flexible  multibody  equations.  In  addition,  we 
describe  an  efficient  Newton-type  iteration  for  dealing  with  the  high-frequency  oscil¬ 
lations  that  arise  from  small  elastic  deformation.  Numerical  experiments  illustrating 
the  effectiveness  of  the  DAE-a  method  are  given  in  Section  4. 

2  Numerical  integration  of  flexible  multibody  sys¬ 
tems 

Based  on  the  spatial  discretization  of  the  elastomechanical  PDE,  the  dynamic  equa¬ 
tions  of  the  response  of  a  discrete  model  of  a  structure  are  given  by 

M^ud-C^u-\-K^u  =  f\t)  (2.1) 
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where  u  €  -K"“  is  the  nodal  displacement,  /®(t)  is  the  time-dependent  function  of 
applied  load,  and  Af®,  C*  and  iif*  are  the  mass,  damping,  and  stiflFness  matrices, 
respectively.  High-frequency  vibrations  usually  result  from  the  spatial  discretiza¬ 
tion.  To  remove  undesired  oscillation,  the  a-method  and  its  variants  are  designed  to 
achieve  controllable  numerical  dissipation  and  A-stability.  These  properties  are  also 
important  in  the  numerical  solution  of  flexible  multibody  systems.  In  what  follows, 
we  formulate  an  extension  of  the  H HT-oc-uiethod  to  solve  the  DAEs  of  constrained 
multibody  systems. 

The  basic  form  of  constrained  multibody  equations  of  motion  is  given  by  (1.1), 
which  is  a  DAE  of  index-3  (see  [1]  for  the  definition  of  index).  Due  to  the  problems 
of  numerical  instability  in  solving  index-3  DAEs,  most  of  the  solution  techniques  for 
(1.1)  have  been  developed  using  differentiation  of  the  constraints  (1.1b).  The  analytic 
solution  of  (1.1)  satisfies 


M{q)q  +  G'^\-  f{q,p,t) 

=  0 

(2.2a) 

G{q)q+  a,  9  +  3(2 

=  0 

(2.2b) 

GiiH+Yt 

=  0 

(2.2c) 

9(9,^) 

=  0 

(2.2d) 

where  (2.2c)  and  (2.2b)  are  the  velocity  and  acceleration  constraints,  respectively. 
They  are  obtained  by  once  and  twice  differentiating  the  position  constraint  (2. 2d). 

Solving  (2.2a)  and  (2.2b)  for  v,  the  underlying  ODE  of  (2.2)  is  obtained 

q=^{q,q,t)  =  M-'{f-G^X)  (2.3) 

where  A  =  f  —  q)  and  q  =  “(^9  +  2^1;  +  fjl)-  For  simplicity 

of  notation,  we  assume  that  M  is  invertible.  Any  conventional  numerical  integration 
method  can  be  applied  to  (2.3),  however  the  numerical  solution  will  not  preserve 
the  constraints  (2.2c)  and  (2.2d).  To  enforce  the  constraints,  the  numerical  solution 
can  be  projected  onto  the  constraint  manifold.  For  an  approximation  q  of  q{t),  q  is 
projected  onto  the  manifold 

M{q){q-q)  +  G'^{q,t)i>  =  0  (2.4a) 

g{q,t)  =  0.  (2.4b) 

This  nonlinear  system  can  be  solved  using  Newton-type  methods,  since  its  Jacobian 
is  nonsingular  near  the  constraint  manifold  [10,  15].  For  the  velocity  constraints,  the 
approximation  v  is  projected  onto  the  velocity  constraint  manifold  by  solving  the 
linear  system 

M{q){v  —  v)  +  G'^{q,t)}j,  =  0  (2.5a) 

G{q,t)v  =  0.  (2.5b) 
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Remark  2.1  The  projection  of  q  onto  the  constraint  manifold  has  been  proposed  by 
many  authors  [18,  23].  A  variant  of  (2.4)  was  used  in  the  half -explicit  extrapolation 
method  [15], 

M{q){q  -q)  +  G^{q,t)i>  =  0, 
where  q  is  obtained  by  extrapolation. 

Applying  the  a-method  (1.3)  to  (2.3),  and  projecting  the  numerical  solution  to 
(2.4)  and  (2.5)  yields 


0„+i  -  (1  +  ci)<f>n+l  +  ot]>n  =  0  (2.6a) 

^{Qn+l)iQn+l  —  qn+l)  +  G^{qn+l)5n+l  =  0  (2.6b) 

Af(?n+l)(Un+l  —  U„+i)  +  G^(9„+i)/I„+1  =  0  (2.6c) 

G{qn+i)vn+i  =  0  (2.6d) 

9{Qn+i)  =  0  (2.6e) 

where 

gn+i  =  9n  +  hvr,  +  /i^((^  -  /3)a„  +  ^SUn+l)  (2.7a) 

Vn+1  =  Un  +  ^((1  -  7)an  +  TOn+l)-  (2.7b) 


For  simplicity,  we  assume  that  the  kinematic  constraints  do  not  depend  explicitly  on 
time,  e.g.,  scleronomic  constraints.  The  accuracy  and  stability  of  (2.6)  can  be  shown 
by  the  standard  analysis  of  (1.3)  and  Gauss-Newton  iteration  of  the  constraints. 

Theorem  2.1  Suppose  g  is  smooth,  e.g.,  g  €  C^,  C?  =  ||  has  full  row  rank  on  the 
constraint  manifold,  and  GMGf^  is  invertible.  Let  c^  be  the  Lipschitz  constant  of  the 
function  <[  of  (2.3),  i.e., 

U{a)-m\<cAW-hl 

There  exists  ho  such  that  for  h  <  ho  the  solutions  ^n+i  owd  of  (2.6)  satisfy 

Ikn+l  -  9(<n+l)ll  =  0{h^),  ||u„+i  -  U(<n+l)||  =  0{h)  (2.8) 

given  consistent  starting  values  {qo,vo,ao,Uo,Po)- 

Proof.  Applying  the  numerical  integration  of  (1.3)  to  the  underlying  ODE  of  (2.2), 
we  denote  =  kn+i,  ^Un+i,  ^^On+i],  where 

On+l  =  (1  +  Q!)<^(9n+l,ti„+i,<n+l)  “  a^(g„,  t„), 

and  qn+i  and  {>„+i  are  obtained  from  (2.7)  using  Cn+i.  Linearization  of  the  HHT 
a-method  yields 

^n+i  =  A„Xn  ,  n  €  {0,l,2,....,iV}.  (2.9) 
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where  Ai  for  i  €  ■{05152,  is  the  loccl  QTnplificQtion  mcitrix.  For  a.  linear 

function  <f>  =  standard  analysis  [11]  shows  that  the  local  truncation  error  of 

X„+i  is  0{h^),  with  a,  7  given  in  (1.3).  Since  is  Lipschitz,  eigenvalue  analysis 
for  the  linear  model  function  can  be  applied  to  the  nonlinear  case  by  the  linearization 
of  (2.9),  see  [2]  for  details. 

Since  GM(^  is  nonsingular  and  G  has  full  row  rank,  locally  there  exists  P{q) 
such  that  P{q)G^{q)  =  0.  Premultiplying  (2.6b,2.6c)  by  Pn+i  =  P(9n+i)  yields 

/  h+lM{qn+l){qn+l  -  qn+l)  ^  =  Q, 

V  Pn+lM{qn+l){Vn+l  “  Un+l)  / 

which  is  a  state-space  ODE  of  (1.1).  The  state-space  solution  together  with  the 
constraints  (2.6d,2.6e)  is  well-defined.  In  addition,  the  local  error  of  {qn+i,hvn+i)  is 
the  same  order  as  that  of  {qn+i,  hvn+i),  provided  that  the  nonlinear  system  is  solved 
with  an  error  tolerance  which  is  less  than  the  truncation  error.  For  more  details 
on  the  proof,  we  refer  to  [24]  (Theorem  3,  pp.  559-561).  Since  Xn+i  is  obtained 
from  the  second-order  a-method,  this  implies  that  the  local  error  of  [^„^i, is 
0(/i^),  provided  that  the  nonlinear  equations  are  solved  with  an  error  tolerance  0{h^). 
Denoting  X„+i  =  [gn+15  hvn+i,  ft^fln+i],  we  conclude  that  the  local  truncation  error  of 
X„+i  is  0{h^).  It  follows  that  (2.8)  holds  for  h  <  ho,  where  ho  depends  on  and  on 
the  Jacobian  and  Hessian  of  g.  □ 

Remark  2.2  The  above  results  hold  for  many  of  the  numerical  integration  methods 
of  structural  dynamics,  including  the  Newmark  ^ -method  [17],  and  a  variety  of  the 
a-methods  [12,  28,  3]. 

Note  that  the  projection  can  be  replaced  by  P„+i  =  P{qn+i)  in  the  above 

analysis  [18].  Therefore,  we  may  use  M{qn+i)  and  G^{qn+-i)  in  (2.6b,2.6c)  instead  of 
M{qn+i)  and  G^(g„+i). 

Remark  2.3  For  the  one-step  a-methods,  the  error  of  the  multipliers  {un+i, Pn+i) 
does  not  depend  on  accumulated  errors  of  the  constraints  from  the  past  state  variables. 
Consider  the  expansion  of  the  constraints  g{qn+i)  about  qn 

fl'(9n+l)  =  fif(9n+l)  +  G{qn+l){qn+l  “  9n+l  +  ^n+l) 

for  some  ||(?n+i^n+ill  =  0(|l9n+i  -  9n+i|P).  Using  the  projection  (2.4),  the  associated 
multipler  Un+i  is  given  by 

K+i  =  {Gn+iM-l,Gl+,)-\giM  -  Gn+iUi)-  (2-10) 

Similarly,  the  projection  (2.5)  yields  the  multiplier 

(2.11) 
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Equations  (2.6)  illustrate  a  straightforward  application  of  the  a-method  to  the 
constrained  multibody  systems.  However,  the  computational  cost  of  this  approach 
is  relatively  large  for  moderate-  to  large-scale  multibody  systems.  Each  iteration  of 
(2.6)  requires  the  inversion  of  the  mass-inertia  matrix  M  and  the  formation  of  /  and 
decomposition  of  GM~^(^  in  the  computation  of  ^n+i,  e-g-,  the  underlying  ODE 

(2.3) .  To  remedy  the  high  cost  of  (2.6),  we  retain  the  Lagrange  multipiler  A  in  (1.1) 
to  avoid  the  inversion  of  large  matrices.  This  is  achieved  because  the  linear  algebra 
involved  in  the  computation  of  <f>n+i  is  the  same  as  those  in  the  projection  steps  of 

(2.4)  and  (2.5).  It  will  be  shown  in  the  next  section  how  this  approach  leads  to  the 
DAE  a-method  which  applies  the  discretization  schemes  of  the  a-methods  in  linear 
structural  dynamics  to  the  index-2  DAE  of  the  constrained  multibody  equations  of 
motion. 


3  The  DAE  a-method 


One  can  reduce  the  computational  cost  of  (2.6)  by  eliminating  the  explicit  solution 
for  ^n+i.  Substituting  the  right-hand  side  of  (2.3)  into  (2.6a)  yields 

On+l  =  (1  +  Q!)Af“+l(/»*+l  ~  ^n+l^n+l)  “  (3.1) 

where  =  M“^(/„  —  G^A„).  Substituting  (3.1)  into  (2.7)  yields 

9n+i  =  9n  +  hv„  +  /i^[(^  -  /?)a„  -  ySa^n]  -I-  (/n+i  -  A„+i  X3.2a) 

Vn+l  =  +  ft[(l  -  7)«n  “  +  /l7M~^l(/„+l  -  G^^.iA„+i)  (3.2b) 

where  ^  a)  and  7  =  7(1  -f  a).  Substitute  (3.2)  into  (2.6b,  2.6c)  to  obtain 

the  DAE  a-method 


Afn+l  (9*1+1  9*1)  ^n+l  d*  f?*i+l  *^*i+l  —  0 

Mn+l{Vn+l  -  K)  -  7^/n+l  +  G^+i/i„+l  =  0 

^n+l^n+1  “  0 

9(9n+l)  =  0 

where 

9t>  =  9n  +  hVn  +  -  fi)an  -  Mn] 

Vn  =  Un  +  ^[(1  -  7)a*i  -  7a<l>n], 


(3.3a) 

(3.3b) 

(3.3c) 

(3.3d) 


(3.4a) 

(3.4b) 


and  Co  =  <^o^and  an  =  (1  -|-a)(^„  4-a^,i_i  for  n  >  1.  Note  that  Un+i  and  g,n+i  in  (3.3) 
comprise  h^pXn+i  and  the  corresponding  correction  terms  by  the  projections. 

Comparing  (3.3)  to  (2.6),  the  solutions  of  the  state  variables  qn+i  and  z;„+i  differ 
by  the  use  of  two  slightly  different  projections.  The  solution  of  (3.3)  is  obtained  from 
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the  projection  along  the  direction  of  while  the  solution  of  (2.6)  is  obtained 

form  the  projection  along  the  direction  of  Mn+i^+i*  On  a  smooth  manifold,  e.g. 
g  e  the  distance  between  the  solutions  can  be  bounded  by  0{h)  for  h  <  ho, 
[18,  24].  However,  the  algebraic  variables  and  /Xn+i  contain  also  the  scaled 
Lagrange  multipier  A  of  (1.1)  which  was  not  accounted  for  by  z>„+i  and  in 

(2.6b)  and  (2.6c),  respectively.  Thus,  convergence  analysis  of  the  algebraic  variables 
(i/„+i,  /i„+i)  remains  for  the  accuracy  and  stability  of  the  DAE  a-method  (3.3).  It  can 
be  shown  via  the  results  for  the  application  of  the  stiffly  stable  numerical  integration 
methods  to  index-2  DAEs  [9]. 

Theorem  3.1  Let  the  assumptions  of  Theorem  2.1  hold,  and  the  initial  value  {qo,vo) 
at  to  satisfy 

11,0  -  «(<„)||  =  0{h^)  ,  11^0  -  »(io)ll  =  0{h^) 

and  g{qo)  =  0,  G{qo)vo  =  0.  The  corresponding  initial  values  of  the  multipliers  are 

I/O  =  ^h  Xq  ,  fiQ  =  7/1A0, 

where  Aq  =  A(to)  of  (1.1).  There  exists  ho  such  that  for  h  <  ho  the  solution  of  (3.3) 

\\qn-q{tn)\\^0{h^)  ,  |K-t;(tn)|l  =  0(/l) 

for  all  n  €  {0, 1, 2, ....,  N}. 

Proof.  The  algebraic  variables  of  (3.3)  can  be  written  as 

t'n+l  =  i^n+l  +  ^h^Xn+l , 

/^n+1  =  An+l  +  7^A„+i 

where  Un+i  and  /in+i  are  due  to  the  projections  of  qn+i  and  Un+i  onto  the  correspond¬ 
ing  constraint  manifolds,  i.e.,  (2.4)  and  (2.5),  respectively.  Using  the  compact  form 
of  (2.9),  e.g.,  =  AnXn,  the  discretization  of  (3.3)  can  be  expressed  by 


9n+l  =  +  ^n+l  +  G„^.iZ/n+l 

(3.5a) 

^  T  ^ 

Wn-l-l  =  Vn  +  han  +  ^^"+1  +  ^n+l^n-|-l 

(3.5b) 

1  - 

On+l  =  On  +  -^On+1 

(3.5c) 

where 

=  ^^^[(1  d"  Oc){M,^^l{fn+l  —  (?^+iAn+l)  —  ~ 

Applying  the  test  equation  q  =  —uj^q  to  (3.5),  three  eigenvalues  of  the  amplification 
matrix  An  as  hu^  —*■  00  are  of  the  magnitudes  [2] 

^  a  1  —  a  1  —  a., 

U  H-a’  1  H-a’  1 
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for  ^  and  7  =  |  —  a.  Thus,  the  a-method  is  stiffly  stable  for  all  a  6 

[—5,0).  Convergence  of  stiffly  stable  integration  methods  for  index-2  DAEs  is  well- 
known  [9].  These  results  yield  the  convergence  of  the  algebraic  variables  in  (3.5).  In 
addition,  for  the  special  case  a  =  0,  i.e.,  the  implicit  trapezoidal  rule,  the  eigenvalues 
corresponding  to  (^n+i,  are  1  and  0.  Convergence  analysis  has  also  been  carried 

out  for  this  case  (for  example  see  [13]). 

Using  these  results  and  Theorem  2.1,  the  local  errors  of  Un+i  and  pn+i  can  be 
obtained  from  (2.10)  and  (2.11),  respectively.  Consequently,  (3.1)  implies  that 


||An+l  -  A(Wi)||  =  0{h) 

where  A(t„4.i)  is  defined  by  (2.3).  For  ||Xn  —  X(<„)||  =  O(h^),  we  conclude  that 

II  -  ,((„„)ll  =  0{h^)  ,  ||8„„  -  .;(t„„)ll  = 

where  qn+i  and  u„+i  are  the  solutions  of  (3.2).  Thus,  the  solution  of  (3.3)  satisfies 
-  9(Wi)||  =  Oih^)  and  ||t;„+x  -  n(in+i)ll  =  Oih),  Vn  €  0, 1, ...,  iV  -  1.  □ 


3.1  Coordinate-split  formulation 

Alternatively,  we  can  directly  apply  the  a-method  to  an  index-1  form  of  (1.1)  using 
nonlinear  projection  operators,  e.g.,  an  annihilation  matrix  P{q)  such  that  P{q)G'^{q)  — 
0.  Premultiplying  (3.3a)  and  (3.3b)  by  such  a  matrix  yields 


P(9„+i)(M„+i(9n+i-9n)-/3/iVn+i)  =  0  (3.6a) 

P(9n+l)(A/n+l(Un+l  -  Wn)  -  Wn+l)  =  0  (3.6b) 

=  0  (3.6c) 

9i<ln+i)  =  0.  (3.6d) 


Accuracy  and  stability  of  the  a-methods  for  ODE  are  preserved,  since  the  algebraic 
variables  have  been  annihilated. 

There  is  a  potential  gain  in  efficiency  for  this  formulation  due  to  the  size- reduction 
of  the  nonlinear  system,  compared  to  (3.3).  An  important  practical  consequence  of 
(3.6)  is  that  {v,p,)  have  been  eliminated  from  the  DAE.  Thus,  the  Newton  iteration 
convergence  test  in  a  numerical  implementation  of  (3.6)  no  longer  needs  to  include 
those  multipliers,  which  can  cause  problems  in  the  direct  numerical  solution  of  (3.3). 
One  could  in  principle  also  consider  removing  {v,p)  from  the  Newton  convergence 
test  in  the  solution  of  (3.3),  however  it  is  not  usually  possible  to  justify  this  action. 
Elimination  of  these  variables  from  the  Newton  convergence  test  in  the  solution  of 
(3.3)  can  lead  to  a  code  which  sometimes  produces  incorrect  solutions.  It  is  the  fact 
that  multiplying  by  the  nonlinear  P(q)  eliminates  (i/,/t)  from  the  nonlinear  system, 
which  allows  these  variables  to  be  excluded  from  the  tests  in  the  solution  of  (3.6). 
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However,  the  computation  of  the  matrix  P{q)  and  its  derivative  can  be  costly  if  not 
impossible.  We  have  dealt  with  this  problem  by  an  efficient  Newton-type  iteration  of 
(3.6)  using  a  family  of  the  C^-matrix  P{q)  [26]. 

To  obtain  a  cheap  representation  of  P{q)  and  its  derivative,  the  CS^-operator  is 
computed  by  the  solution  of  a  class  of  pseudo-inverses  of  the  constraint  Jacobian 

G(,). 

Definition  3.1  [Coordinate-Splitting  Matrix]  Let  X  and  Y  be  the  matrices  whose 
columns  constitute  the  standard  basis  such  that  ||(G(g)y)“^  ||  is  bounded  in  a  neigh¬ 
borhood  Uo  of  qo.  The  coordinate-splitting  matrix  is  defined  by 

Piq)  =  -  QiqfY-^  =  X^il  -  G(qnG{q)Y)-'^Y^)  (3.7) 

where  Q(q)  =  {G(q)Y)-'a{q)X . 

Remark  3.1  From  the  construction  of  the  CS  matrix  P{q),  one  can  easily  see  that 
P{q)G^{q)  =  0  for  all  q,  i.e.,  P{q)  is  orthogonal  to  range{G^).  Furthermore,  the  row 
vectors  of  P{q)  are  orthonormal,  i.e.,  P(q)^P{q)  is  the  identity. 

The  computation  of  P{q)  can  be  carried  out  using  the  LC/-factorization  of  the 
constraint  Jacobian  matrix.  Applying  Gaussian  elimination  with  row-pivoting  to  G^ 
yields 

Err,---ErG^  =  Lrr,---L^U  (3.8) 

where  E,-  is  the  elementary  permutation  and  Li  is  a  Gauss  transformation,  i  €. 
{1, 2, ...,  m}.  From  the  factorization,  we  have 

[Y,X]  =  E  =  E,n--Ei.  (3.9) 

Using  the  standard  solution  technique  by  TU-decomposition,  the  projected  vector 
P{q)r  can  computed  in  a  straightforward  and  relatively  cheap  way.  In  addition,  the 

derivative  can  be  computed  by  the  same  factorization  of  G^  and  the  interme¬ 

diate  result  s  =  —{GY)~^Y^r  from  the  computation  of  P{q)r. 

Remark  3.2  Alternatively,  one  can  apply  Q R-factorization  to  G^  for  the  computa¬ 
tion  of  P{q)r.  Using  Q R-factorization  with  partial  column  pivoting  [7],  we  obtain 

=  (3.10) 

where  Ei  is  the  elementary  permutation,  Li  is  the  Householder  matrix,  i  €  {1,2,  ...,m}, 
and  U  is  upper  triangular.  The  last  (n  —  m)  columns  of  the  orthogonal  matrix 
L  =  Li  constitute  a  basis  for  the  null-space  of  G.  Thus  we  can  write 

y,^]  =Z  =  Xl•••Z„^.  (3.11) 

Note  that  X  and  Y  are  usually  subsets  of  the  standard  basis  in  EC . 
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An  efficient  implementation  of  Newton-type  iteration  for  (3.6)  has  been  presented 
in  [26],  where  the  Jacobian  of  the  nonlinear  system  was  derived.  More  importantly, 
a  modification  of  the  Jacobian  was  obtained  to  enhance  the  convergence  of  Newton- 
type  iterations  for  high-frequency  oscillations.  This  is  achieved  by  eliminating  the 
term 

dq  dq 

where  s  =  —{GY)~^Y^r,  from  the  Jacobian.  This  term  represents  the  sensitivity 
of  the  potential  force  due  to  the  constraints.  When  the  solution  is  not  close  to  an 
equilibrium  on  the  constraint  manifold,  this  term  can  be  dominant  in  the  Jacobian 
of  the  nonlinear  system.  Consequently,  the  iterative  solution  is  forced  to  resolve 
the  high-frequency  oscillations  if  the  constraint  violation  is  large  compared  to  the 
amplitude  of  the  high  modes.  In  highly  oscillatory  flexible  multibody  systems,  this 
problem  can  occur  at  larger  time  integration  stepsizes.  Using  the  modification  of  the 
CS  iteration,  abbreviated  CM,  the  DAE  a-method  in  the  form  of  (3.6)  permits  a  fast 
convergence  of  the  solution.  In  addition,  the  computational  cost  of  the  C'5'-matrix 
is  insignificant  compared  to  one  iteration  of  (3.3)  for  DAEs  with  a  small  number  of 
constraints,  such  as  large-scale  flexible  multibody  systems.  We  will  present  several 
numerical  examples  involving  high-frequency  vibrations  to  illustrate  the  advantages 
of  the  CM  iteration  in  the  next  section. 


4  Numerical  examples 


The  DAE  HHT  a-method  in  (3.3)  can  be  used  for  any  numerical  integration  scheme 
for  linear  structural  dynamics.  For  the  numerical  solution  of  the  test  problems,  we 
will  use  the  DAE  formulas  which  are  extensions  of  the  generalized  a-method  [3].  The 
generalized  a-method  for  (1.2)  is  given  by 


(1  ®m)®n+l  "t"  —  (1  ^}fni 


(4.1) 


where  the  position  and  velocity  states  are  defined  by  (1.3b)  and  (1.3c),  respectively. 
The  parameters  of  the  method  are  obtained  by  a  user-specified  value  of  high-frequency 
dissipation  p  €  [0, 1]  (i.e.,  p  =  0  maximum  dissipation,  p  =  1  none),  such  that 


Olm  — 


2^-1 
/>  +  ! 


aj 


P 

/>+!’ 


1 

7  =  2  “  +  a/  , 


^  ~ 

^  4 


This  yields  a  second-order,  A-stable  family  of  numerical  integration  methods.  Using 
(4.1)  in  (3.3)  or  (3.6),  we  need  only  replace  (3.4)  by 


9n  =  4n  +  +  t’|(i  -  — 

L  1  OLffi  1  OLm. 


(4.2a) 
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(4.2b) 


Vn  =  fn  +  ^[(1  -  ~  I  — a 

i  —  OLm  -L 


and  ^  =  /?f^’  generalized  a-method  is  parametrized  by  the  size 

of  the  dominant  eigenvalue(s)  p  of  its  amplification  matrix  A{—{hu}Y)  as  hu)  —*  0.  It 
gives  the  most  comprehensive  range  of  the  user-specified  value  of  the  spectral  radius 
in  the  high-frequency  limit,  while  low-frequency  dissipation  is  minimized  [3]. 


4.1  Stiff  spring  pendulum 

To  demonstrate  the  user-specified  high-frequency  dissipation  in  the  generalized  a- 
method,  we  use  the  stiff  pendulum  example.  In  Cartesian  coordinates,  a  simple  stiff 
pendulum  model,  with  unit  mass  and  gravity,  may  be  expressed  as 


0 

=  X  —  u 

(4.3a) 

0 

=  y-v 

(4.3b) 

0 

=  it  -f  xA 

(4.3c) 

0 

=  i)  +  y\  —  1.0 

(4.3d) 

e^A 

+  y'^  —  1.0 
y/x^  -(-  y'^ 

(4.3e) 

where  the  stiff  spring  of  natural  length  1.0  and  stiffness  ^  is  attached  to  the  center 
of  mass  of  the  pendulum.  Denote  by  h  the  stepsize,  6  the  maximum  deflection  of  the 
spring,  and  the  potential  energy  of  the  system 


Then  the  frequency  u;  of  the  system  is  proportional  to  l|dy||  such  that 


1  r±_ 

e\ll+6 


<  a;  <  - 
c 


As  e  — >  0,  the  dominant  pair  of  eigenvalues  approaches  ±oo  along  the  imaginary  axis. 
This  problem  is  discussed  in  more  detail  in  [25]. 

The  numerical  solution  of  the  index-1  DAE  (4.3)  has  been  carried  out  applying  the 
generalized  a-method  to  (3.3).  Using  the  initial  values  (xo,  j/o,  wo,  t^o)  =  (0.9, 0,0,0), 
e2  =  10-3,  and  Newton  iteration  with  tolerance  10"^°,  the  results  of  several  combi¬ 
nations  of  the  high  frequency  dissipation  p  and  steplength  h  for  a  0  to  5  second  sim¬ 
ulation  are  given.  The  total  energy  is  plotted  in  Figure  4.1,  with  the  high-frequency 
dissipation  parameter  p  varying  from  1  to  0.  The  solution  trajectory  on  the  coordi¬ 
nate  plane  is  shown  in  Figure  4.2  with  the  stepsize  varying  from  10'^  to  10”^  The 
results  illustrate  the  effective  damping  of  the  DAE-a  method.  As  p  — »  0  or  the  step- 
size  h  is  increased,  the  numerical  solution  approaches  the  equilibrium  of  the  spring 
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while  maintaining  the  slow  swing,  e.g.,  see  Figure  4.2.  Since  the  accuracy  of  the  nu¬ 
merical  solution  is  second-order,  the  increasing  error  in  the  total  energy  corresponds 
to  increasing  h  in  the  picture.  The  controllable  damping  of  the  DAE  a-method  is 
demonstrated  by  the  stiff  spring  pendulum  example.  However,  the  stepsize  may  be 
restricted  to  the  convergence  of  the  nonlinear  equations  resulting  in  the  discretization. 
In  the  next  subsection,  we  present  a  nonlinear  bushing  example,  which  illustrates  that 
the  CM  iteration  of  (3.6)  achieves  better  convergence  than  solution  of  (3.3)  via  New¬ 
ton  iteration,  thus  allowing  effective  damping  of  high-frequency  oscillations  by  taking 
larger  time  steps. 


Figure  4.1:  Stiff  Spring  Pendulum  Example;  =  10  TOL  =  10 


Trajectory  on  X-Y  Plane:  rho^.O,  h>0.00 1,0.0 1,0.1 


Figure  4.2:  Stiff  Spring  Pendulum  Example;  =  10  TOL  =  10 
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4.2  Nonlinear  bushing  force 

The  next  example  is  a  nonlinear  three- variable  bushing  example.  Using  2-D  Cartesian 


coordinates,  the  problem  may  be  represented  as 

0  =  X  —  —fx  (4.4a) 

0  =  y-^/y-1  (4-4b) 

0  =  0  + j^^-^(cos0(^/v-M)-sin0^/r)  (4.4c) 

0  =  x^  +  y^-1  (4.4d) 


where  \  - x +  \cose,  fy  =  -y sin $,  and  e  =  10-^  This  problem  was  solved 
using  both  the  stabilized  index-2  DAE  (3.3)  and  the  CM  iteration  of  (3.6).  Both  the 
formulations  work  well  for  the  problem.  The  CM  iteration  permits  a  larger  stepsize 
and  therefore  a  much  quicker  dissipation.  Some  results  are  compared  in  Figures  4.3 
and  4.4.  With  the  initial  values  x  =  0.8  ,  y  =  0.6  and  6  =  0.0  and  the  tolerance 
for  the  Newton  convergence  10“^,  the  figures  show  the  response  of  the  9  variable  for 
different  parametric  values.  The  simulation  was  run  for  T  =  407r  x  10  ®  seconds,  i.e., 
roughly  6  cycles  of  9  and  around  20  cycles  for  x  and  y. 

ButNng  Fore*:  h  m  T500,  rtw  ■  0.1 


Figure  4.3:  Bushing  Example,  Stabilized  Index-2  Formulation;  e  =  10  TOL  —  10  ^ 

( 1  ^  /  ic  0  \  ^ 

The  CS  iteration  is  implemented  using  P  =  (  ^  ^  1  /  ’  x  thI  during 

the  simulation.  In  Table  4.1,  we  list  the  damping  measurement  of  the  solution  of  the 
CM  iteration  and  that  of  (3.3).  It  is  measured  by  the  ratio  of  the  initial  energy  and  the 
energy  retained  at  the  end  of  the  simulation.  Because  a  straightforward  measurement 
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Bu«Hng  Forot :  h  »  T/500,  rtK>  -  0.1 


Figure  4.4:  Bushing  Example,  CM  Iteration  for  500  Timesteps;  e  =  10  TOL  =  10  ^ 


No.  steps 

p  =  0  (CM) 

p  =  o 

mstmi 

WOi 

P  =  i  (CM) 

■muM 

■IlillM 

■WiSM 

mmum 

■IIWM 

■IIBIM 

max.  ho 

T 

_ 40 _ 

mm 

T 

_ 20 _ 

mm 

T 

_ 80 _ 

wmm 

Table  4.1:  Final  Energy  vs.  Initial  Energy  for  the  Bushing  Problem;  T  =  407r  x  10“® 

of  damping  from  the  displacement  response  is  difficult,  the  energy  norm  has  been 
conveniently  used.  The  last  row  of  Table  4.1  shows  the  maximum  stepsize  that  can 
be  taken  without  causing  a  failure  in  Newton  convergence.  For  this  example,  the 
CM  iteration  allows  stepsizes  twice  as  large  as  those  of  the  index-2  formulation.  In 
addition,  the  CM  solution  exhibits  more  rapid  convergence  towards  the  slow  solution, 
e.g.,  approaching  a  smaller  energy  level  than  that  of  (3.3).  This  is  illustrated  in  Table 
4.1,  where  the  same  h  and  p  were  used  for  both  CM  and  (3.3)  solutions. 

4.3  Flexible  slider-crank 

A  flexible  slider-crank  mechanism  is  presented  schematically  in  Figure  4.5.  The  mech¬ 
anism  consists  of  a  rigid  crank  connected  with  a  flexible  rod  driving  a  slider  on  a 
straight  line.  Three  coordinates  (<f>i,<l>2,X3)  are  used  for  the  rigid  body  equations 
of  motion,  and  the  flexible  rod  is  modeled  by  a  finite  element  approximation  of  a 
linearized  Bernoulli-Euler  beam.  For  the  details  of  this  model,  we  refer  to  [22]. 

Using  a  simple  beam  obtained  from  the  linearization  with  respect  to  the  vertical 
displacement  by  setting  the  longitudinal  displacement  to  zero,  the  two  constraint 
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equations  are 


li  sin  <f>i  +  I2  sin  <f>2  +  92fc+i  cos  <^2 
X3  —  li  cos  (f>i  —  I2  cos  <l>2  +  92fc+i  sin  4>2 


) 


where  /i  is  the  length  of  the  crank,  l2  is  the  length  of  the  rod,  qj  are  the  nodal  coor¬ 
dinates  induced  by  the  k  elements  and  92^+1  is  the  nodal  coordinate  at  the  boundary. 
Since  the  constraint  involves  only  the  boundary  nodes,  the  computation  of  P{q)  for 
the  CS  formulation  can  be  reduced  to 


COS  ^2  [  h^OS<f>i 

sin  ^2  )  V  hsm^i 


I2  cos  <f)2  —  92fc+i  sin  <^2  ^  yT 

I2  sin  ^2  +  92Jfe+i  cos  (j>2  ) 


(4.5) 


where  consists  of  the  column  vectors  63  and  62^+1  and  consists  of  ei  and  62 
such  that  tj  is  the  ith  standard  base  of  The  matrix  P{q)  is  applied  directly 

to  the  4x1  vectors  of  Xj{M{q)  —  /)  with  Xj  =  (61,62,63,62^+1).  The  rest  of  the 
2k  +  \  equations  remain  unchanged.  The  computational  cost  of  the  C S'-matrix  P  is 
equivalent  to  that  of  P.  Increasing  the  number  of  elements,  the  computation  of  P 
remains  unchanged. 

A  one-second  simulation  was  run  with  a  driving  torque  of  0.2  N-m  from  0-0.3 
second.  The  initial  velocity  and  all  position  values  are  zero  except  X3  =  0.45.  The 
Newton  iteration  tolerance  was  set  to  10~®.  The  spatial  discretization  of  the  flexible 
rod  was  carried  out  using  k  =  A,  where  the  magnitude  of  the  maximum  eigenvalue 
is  about  1.3  x  10^.  Using  (4.5)  in  the  CM  iteration,  Newton  convergence  can  be 
attained  with  larger  stepsizes  compared  to  the  index-2  DAE  form,  see  Table  4.2.  The 
high-frequency  oscillation  of  the  nodal  coordinates  is  shown  in  Figure  4.6,  where  the 
amplitude  of  the  vibrations  is  small  compared  to  the  rigid  motion.  The  numerical 
solution  corresponding  to  the  slow  motion,  i.e.,  (^i,<^2,®3),  is  preserved  by  the  iter¬ 
ations  of  CM  and  stabilized  index-2  DAE,  as  shown  by  the  (f>2  plots  in  Figure  4.7. 
The  stepsize  h  =  ^  is  about  one  period  of  the  highest  frequency.  Applying  stronger 
damping  to  the  high  frequency  modes,  the  solution  at  the  midpoint  of  the  flexible 
rod  is  presented  in  Figure  4.8  with  p  =  0.5, 0.0  and  h  =  The  high  modes  are 
eliminated  in  the  numerical  solution,  comparing  to  Figure  4.6.  For  various  combi¬ 
nations  of  p  and  h,  the  total  energy  plots  are  presented  in  Figures  4.9  and  4.10.  At 
various  levels  of  damping  corresponding  to  the  different  values  of  p  and  h,  only  the 
higher  modes  are  removed  in  the  numerical  solution.  It  is  shown  in  Figure  4.11  that 
faster  numerical  dissipation  of  the  high  modes  is  achieved  via  the  CM  iteration,  where 

h  =  =  0-75. 
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Table  4.2:  Maximum  Stepsize:  CM  vs.  Index-2  DAE,  p  =  0.75 


Longitudinal  Displacement 


Figure  4.5:  Slider-Crank  Mechanism  with  a  Flexible  Rod 
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Figure  4.6:  Nonlinear  Oscillations  in  Nodal  Coordinates 


nwdbtoSMarCrwtk  Example  h- 1/500,  rho-0.75 


Figure  4.7:  Rigid  Motion  of  Coordinate  <f>2 
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Figure  4.8:  Damped  Vibration  on  Midpoint  of  Flexible  Rx)d 


SNdtr-Crank  Examplt:  h-1AOO.  rho  -  0.0.5.075.1 


Figure  4.9:  Total  Energy  Comparison,  p  =  0,0.5,0.75,1.0 
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S«d*r-Crtr*  Exvnpk;  h  - 1/300, 1/500, 1/750,  rtio-0.5 


Ttm* 


Figure  4.10:  Total  Energy  Comparison,  ^  =  3555  sm’  tIo 


ntiMte  8lktor-Gr«*  Eumfsta:  h  •  1/500,  rho-0.7S 


Figure  4.11:  Total  Energy  Comparison,  CM  vs.  Stabilized  Index-2  DAE 
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